Informe Hito 3: Trabajo con datos sismológicos

Introducción

Los sismos en Chile son fenómenos muy recurrentes e impredecibles. El cinturon de fuego del pacifico hace que paises que contiene tengan gran frequencia de sismos. A diferencia de decadas anteriores el numero de sismos detectados incremento significativamente gracias a nuevas estaciones sismico-detectoras instaladas en Chile.

Los sismos puden tener consecuencias de alto costo humano y economico, por ello es importante no solo intentar predecir su ocurrencia sino también estudiar el efecto que podrían provocar el sismo en otros ámbitos. Es decir el impacto que tiene en otras areas. En este informe nos centraremos en los recursos humanos, por el accesso a estos datos que permite el Departamento de Estadisticas e Información de Salud (DEIS) en Chile.

Se investiga si eventos sismicos tiene relación con las Enfermedades Transmitidas por Alimentos (ETA) como tambien impacto al incremento de Ingresos a Hospitales (IH). El poder determinar si existe una relación entre estos eventos podría ayudar a desplegar infraestructura y recursos humanos después de un desastre, además de ayudar en la generación de políticas públicas que ayuden a mitigar los efectos secundarios.

Nota. Desde el Hito 2 se trabajan las hipótesis recién mencionadas, puesto que en el Hito 1 se barajó la posibilidad de basarse en el estudio de la relación entre eventos fluviales y los sismos. Este informe incluye este análisis con fines informativos. Los análisis finalmente indicaron que las relaciones no son claras para los sismos, por lo menos para Chile.

Hipótesis

Hito 1

Se trabajó con dos estudios principales: el primero, sobre la relación de los sismos con las lluvias monzónicas en el Himalaya, y el segundo, sobre la influencia del sol en ellos. Los resultados del primer estudio mostraron que la frecuencia de sismos en el Himalaya se veía disminuida con los eventos monzónicos. Del segundo estudio, se extrajo que los datos sugerían que la mayor cantidad de sismos ocurrían en horarios diurnos.

La expoloración de los datos sugirió que no habían relaciones claras entre los eventos estudiados como para generar un predictor. En este entendido, se añade al conjunto de dataset una base de datos de egresos hospitalarios, con el fin de direccionar la experimentación a hacia el aumento o la disminución de estos, incluyendo de la misma forma las enfermedades transmitidas por alimentos.

Hito 2

Se trabajó con los datos de enfermedades transmitidas por alimentos e ingresos hospitalarios. Se planteraon asimismo hipótesis todavía vigentes:

Configuración A: Correlación entre sismos y cantidad de egresos hospitalarios en Chile.

Hipótesis: Los sismos inciden en la cantidad de ingresos hospitalarios por cada región. Se llegó a que efectivamente existía cierta correlación, pero faltaba indagar más.

Configuración B: Correlación entre sismos y cantidad de reportes de enfermedades transmitidas por alimentos en Chile.

Hipótesis: Los sismos inciden en la cantidad de enfermedades transmitidas por alimentos (ETAs), particularmente por deshidratación (como se sugirió con los datos del Hito 1). Se llegó a que en realidad, había una correlación con la cantidad de ETAs reportadas y los sismos, no en particular por deshidratación.

Hito 3

Al tratarse de una extensión y mejora del hito 2, no se realiza cambios en las hipótesis ni configuraciones realizadas en el Hito 2.

Metodología

Se busca juntar los datasets para buscar las correlaciones mencionadas en las hipótesis. La forma de juntarlos es mediante la posición de los eventos (latitud/longitud) o región, en caso que corresponda, y la fecha en que se produjeron. De esta manera se buscará la correlación entre los eventos planteados para confirmar o rechazar las hipótesis.

Se procedió a continuar con lo realizado en el hito anterior, haciendo las proyecciones especificadas ahí y también nuevas ideas pensadas posteriormente.

Descripción de los datos

Fuente de Datos

Conjunto de dataset actualmente utilizado:

Dataset	Objetivo del dataset
Datos sismológicos	Extraer los datos de sismos y poder compararlos con los egresos y las ETAs
Egresos hospitalarios	Dado un sismo, experimentar el aumento o la disminución de egresos hospitalarios
ETAs	Dado un sismo, experimentar el aumento o la disminución de enfermedades transmitidas por alimentos.

Las tablas de Egresos hospitalarios y de ETAs tienen mayor cardinalidad que la de datos sismológicos.

Datos sismológicos

IRIS Incorporated Research Institutions for Seismology o por sus siglas en español, Instituciones incorporadas de investigación sismológica, provee equipamiento y acceso sísmico y otros datos alrededor del mundo, cortesía de una red sismógrafos de la comunidad científica internacional y de Estados Unidos. Siendo Chile uno de los paises con más frecuencia de sismos en el mundo, IRIS clasifica los sismos dentro de un cuadro limitado por coordenadas definidas por el mismo instituto (máxima latitud=-15.400, mínima latitud=-57.000, máxima longitud=-63.800, mínima lon=-83.100). Como se muestra en el gráfico a continuación:

Ultimos 100 sismos en la zona al rededor de Chile

Dentro de los datos que se utilizaran para el análisis están: longitud, latitud, magnitud, profundidad, y tiempo. Longitud y latitud son parte del sistema de coordenadas geográficas es un sistema que referencia cualquier punto de la superficie terrestre y que utiliza para ello dos coordenadas angulares, latitud (norte o sur) y longitud (este u oeste). Magnitud es una medida que tiene relación con la cantidad de energía liberada en forma de ondas. Profundad define la distancia entre el epicentro de un terremoto con respecto al nivel del mar. Tiempo es el registro del evento.

El gráfico de abajo muestra la cantidad de reportes recolectados por año, dentro de las inferencias, y consultas hechas los outliers que puedan ser registrados con equipos antiguos están exentos por ser mínima en su cantidad.

# load libraries
library(tidyverse)
library(lubridate)

ggplot(season_data, mapping = aes(x=year)) +
  geom_bar()

Cantidad de sismos por año

Datos de registros de ingresos hospitalarios

Los datos de Ingresos Hospitalarios (IH) son bastos, la collección que se pudo recabar data de los años 2018 a 2019. Primero vemos la cantidad de elementos por año es basta y varia, como se ver el grafico de abajo. No encontramos ninguna justificación por la que esto ocurre en la pagina del DEIS.

f <- data.frame(year= c(2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018),
data=c(nrow(egr_2008), nrow(egr_2009), nrow(egr_2010), nrow(egr_2011), nrow(egr_2012), nrow(egr_2013), nrow(egr_2014), nrow(egr_2015), nrow(egr_2016), nrow(egr_2017), nrow(egr_2018)))
plot(f, type="l")

Esta información cuenta con varios campos, dentro de los mas interesantes son los mensionados a continuacion: COMUNA, FECHA_EGR, DIAS_ESTAD, REGION, y DIAG1. Pronto vemos que el numero de diagnosticos unicos son miles por lo tanto no sirven de ayuda para filtrar estos datos.

Dentro de este análisis buscamos encontrar los ingresos hospitalarios. Sin embargo estos los datos originales no muestran este datos, sin embargo podemos genrar una nueva columna de fecha de ingreso utilizando la formula FECHA_EGR - DIAS_ESTAD. Un problema encontrado para generar esta nueva columna fue que los formatos de las fechas puede variar de acuerdo el año que se encuentre. El siguiente codigo muestra el algoritmo que se utilizo para este proceso.

ddd <- data.frame(a=c(1,2,3))
ddd <-mutate(ddd, b= a + 1)

a <- c(1,2,3) + c(1,2,3)
data <- egr_2018

formats <- c("%d-%m-%Y","%d-%m-%Y","%d-%m-%Y","%d-%m-%Y",
            "%Y-%m-%d","%Y-%m-%d","%Y-%m-%d","%Y-%m-%d","%Y-%m-%d","%Y-%m-%d","%Y-%m-%d")

dclist <- list()
counter = 1
for (data in dlist) {
  
  lastDate <-  as.Date("1900-01-01")
  dtime <- as.Date(data[['FECHA_EGR']], formats[counter])
  dtime <-  dtime - data[['DIAS_ESTAD']]
  season_data <- mutate(data,
    year_ing = year(dtime),
    month_ing = month(dtime),
    mday_ing = mday(dtime),
    number_days = as.numeric(dtime - lastDate)
  )
  dclist[[counter]] <- season_data
  counter = counter + 1
}

Datos de enfermedades transmitidas por alimentos

Los datos de salud alimentaria fueron extraídos de la página web del Departamento de Estadísticas e Información de Salud, en www.deis.cl. Se extrajeron datos de los brotes de enfermedades transmitidas por alimentos desde el año 2011 hasta el 2017, siendo estos todos los disponibles.

Para hacer una primera exploración de los resultados, se contó la cantidad de casos de enfermedades transmitidas por alimentos reportadas cada año para algunas regiones afectadas por terremotos grandes que ocurrieron en el periodo 2011-2017, aquellos de Iquique el año 2014 y de Illapel el año 2015.

Ahora, los sismos vienen con las columnas Year, Month y Day, que fueron generados con una función similar a agregarFechas que se encuentra en el anexo. También incluyen las regiones: estas fueron generadas mediante la función cargar_sismos_por_region que se encuentra en el anexo.

Análisis hitos anteriores

Se cargan los sismos con:

sismos <- read.delim("http://anakena.dcc.uchile.cl/~rllull/CC5206/sismos.csv",
                     header = TRUE, sep = ",", quote = "\"", dec = ".", fill = TRUE, comment.char = "#")

El código para las ETAs se simplificó mucho, logrando unir los 7 datasets de ETAs en 1. De esta manera, se cargan los datos de la forma:

# Carga de las ETAs
etas <- read.delim("http://anakena.dcc.uchile.cl/~rllull/CC5206/etas2011_2017.csv",
                   header = TRUE, sep = ",", quote = "\"", dec = ".", fill = TRUE, comment.char = "#")

# Agregarle las fechas a las ETAs
etas <- agregarFechasEtas(etas)

Después, se guardan los datos en data frames cant_etas y cant_etas_deshidratacion, para obtener algunos gráficos:

Reporte de Region por año

Diseño experimental

Dentro del disenio experimental se inclino a realizar un análisis de clustering y clasificación de sismos con respecto a eventos ETA y IH. Este necesito de un riguroso análisis de los datos para poder procesar y preparar la data. Dentro de los cambios mas fuertes fue el de la categorización de los datos sismos (datos geograficos) en regiones a Chile. Ademas del filtrado de los caracteristicas utiles para este análisis.

Se plantea una clasificación de sismos binaria y otra de cuatro dimensiones de incrementos en ETA y IH. En cuanto a los clusterings se muestra la correlación los sismos tienen con los incrementos de ETA y IH, asi encontrar una segementación que nos ayude la relación de estos fenomenos.

:TODO brief conclusion

Limpieza de los datos

Datos sismológicos

Los datos sismicos son la principal fuente de información para los experimentos. Este requirio de los siguiente pasos.

• Con el nuevo objetivo definido, se buscó encontrar la correlación que se produce entre los sismos los informes de ETA e IH. Para ello, se decidió utilizar solamente aquellos sismos de magnitud mayor o igual a 5 Mw (magnitud en la escala de Ritchter).

• Se realizo la clasificación regional de los datos filtrando sismos que pertenezcan a Chile, agregando una columna para agregar el nombre de la region correspondiente (el codigo se encuentra en anexos).

• Se agrego datos en representación de las fechas para facilitar las consultas a la data.

• Por ultimo se filtro aquellas columnas que si representaban un verdadero aporte para los casos de estudio (‘Latitude’, ‘Longitude’, ‘Depth’, ‘COMUNA’).

Nota. Lo anterior mensionado es suficiente para poder utilizar los datos, sin embargo estos no son suficientes como veremos mas adelante. Por ejemplo, es necesario escalar los datos para el clustering y clasificación . Ademas los datos satelite (ETA e IH) obligan que existan ciertas restricciones en cada análisis.

Aumento registrado según la región

Procedemos a hacer un aggregate con tal de contar cuánto aumentó, en promedio, por cada magnitud y por cada región. Además, se puede extraer si se aprecia un aumento según cada magnitud y cada intervalo, obteniendo un resultado en cantidad de ETAs registradas posteriores o anteriores al sismo en promedio.

mean <- function(a) {sum(a)/length(a)}

# Magnitud 5 o más
(agg_5_7d <- aggregate(Aumento ~ Region, mag5_7d, FUN=mean))

# Magnitud 6 o más
(agg_6_7d <- aggregate(Aumento ~ Region, mag6_7d, FUN=mean))

# Magnitud 7 o más
(agg_7_7d <- aggregate(Aumento ~ Region, mag7_7d, FUN=mean))

Se aprecia que en general en promedio se registran aumentos en los reportes de ETAs para cada región. En la parte de Minería de datos se estudia con mayor detalle este hallazgo.

Procesamiento de datos para clasificación (agregacion)

Para generar correctamente un modelo de clasificación, se necesitó pensar en un procedimiento efectivo que dijera con certeza si un terremoto realmente provocó daños. Esto, con base en la hipótesis de que mientras más daños pueda generar un sismo, más informes por ETAs se generan. Además, se debió agregar la dimensión de temporalidad, necesaria para apreciar los datos.

De esta manera, se pensó en lo siguiente: si en los n días antes del sismo se registraron menos informes por ETA que en los n días después, quiere decir que el sismo aumentó los informes y por ende los daños. De lo contrario, los disminuyó. Este valor se denominó Δ, con \(\Delta = etas_{post} - etas_{ant}\).

En el Hito 2 se representó con un 1 si el Δ era positivo, y un 0 si era negativo. El nuevo Δ permitió aumentar la cantidad de clases a estudiar.

Imagen de cómo se calculó aumento

Se extraen para diferentes magnitudes (5, 6 y 7) e intervalos (de 3, 7, 14 y 30 días). Se crea una función para trabajar los datos más fácil.

generar_datos <- function(sismos, etas, magnitud_min, dias, escribir=TRUE) {
  sismos_mag <- sismos[sismos$Magnitude >= magnitud_min
                       & sismos$Year >= 2011 & sismos$Year <= 2017, ]
  # Generar los dataframes
  # datos es para trabajarlo en el árbol de decisión
  datos <- data.frame(Magnitud = rep(0, nrow(sismos_mag)),
                      Profundidad = rep(0, nrow(sismos_mag)),
                      Region = rep(0, nrow(sismos_mag)),
                      Aumento = rep(0, nrow(sismos_mag)))
  # datos_con_fecha es para trabajarlo en lo que viene
  datos_con_fecha <- data.frame(datos,
                              Año = rep(0, nrow(sismos_mag)),
                              Mes = rep(0, nrow(sismos_mag)))
  for (i in 1:nrow(sismos_mag)) {
    sismo <- sismos_mag[i, ]
    # Tomar los dias días anteriores
    etas_ant <- etas[sismo$num_date - dias/30/12 <= etas$num_date
                     & etas$num_date <= sismo$num_date, ]
    # Tomar los 7 días posteriores
    etas_post <- etas[sismo$num_date < etas$num_date
                      & etas$num_date <= sismo$num_date + dias/30/12, ]
    # La diferencia es la cantidad de columnas de las posteriores menos las anteriores
    dif <- nrow(etas_post) - nrow(etas_ant)
    
    # Asignar los valores a los datos
    datos[i, ]$Magnitud <- sismo$Magnitude
    datos[i, ]$Profundidad <- sismo$Depth
    datos[i, ]$Region <- sismo$Region_Number
    # A diferencia del Hito 2, se guarda el Δ completo
    datos[i, ]$Aumento <- dif
    # Guardar datos a retornar
    datos_con_fecha[i, ] <- datos[i, ]
    datos_con_fecha[i, ]$Año <- sismo$Year
    datos_con_fecha[i, ]$Mes <- sismo$Month
    
  }
  # Escribir el csv con los datos obtenidos
  if (escribir) {
    nombre_archivo <- paste("datos_mag", magnitud_min, "_", dias, "d.csv", sep="", collapse=NULL)
    write.csv(datos, nombre_archivo, row.names=FALSE)
  }
  (datos_con_fecha)
}

Para generar los datos, se procedió de la siguiente forma:

# Sismos de magnitud mayor o igual a 5
mag5_3d <- generar_datos(sismos, etas, 5, 3)
mag5_7d <- generar_datos(sismos, etas, 5, 7)
mag5_14d <- generar_datos(sismos, etas, 5, 14)
mag5_30d <- generar_datos(sismos, etas, 5, 30)

# Sismos de magnitud mayor o igual a 6
mag6_3d <- generar_datos(sismos, etas, 6, 3)
mag6_7d <- generar_datos(sismos, etas, 6, 7)
mag6_14d <- generar_datos(sismos, etas, 6, 14)
mag6_30d <- generar_datos(sismos, etas, 6, 30)

# Sismos de magnitud mayor o igual a 7
mag7_3d <- generar_datos(sismos, etas, 7, 3)
mag7_7d <- generar_datos(sismos, etas, 7, 7)
mag7_14d <- generar_datos(sismos, etas, 7, 14)
mag7_30d <- generar_datos(sismos, etas, 7, 30)

Clustering

Clustering ETA

Análisis del aumento con boxplot

Se puede estudiar el balance de los datos por cada magnitud según la cantidad de días de intervalo tomados. Del resultado anterior se sugiere un aumento en la cantidad de ETAs registradas según la cantidad de días de intervalo considerados para los sismos.

par(mfrow=c(1, 3))

outl_mag5 <- boxplot(mag5_3d$Aumento, mag5_7d$Aumento, mag5_14d$Aumento, mag5_30d$Aumento, main="Aumento para magnitud >= 5\ny 3, 7, 14 y 30 días de intervalo")$out
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

outl_mag6 <- boxplot(mag6_3d$Aumento, mag6_7d$Aumento, mag6_14d$Aumento, mag6_30d$Aumento, main="Aumento para magnitud >= 6\ny 3, 7, 14 y 30 días de intervalo")$out
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

outl_mag7 <- boxplot(mag7_3d$Aumento, mag7_7d$Aumento, mag7_14d$Aumento, mag7_30d$Aumento, main="Aumento para magnitud >= 7\ny 3, 7, 14 y 30 días de intervalo")$out
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

Aumento ETA con respecto a Sismos

Se ve que los datos para 30 días tienen mucha dispersión, y aquellos para 3 días están muy centrados al 0. Se aprecia además que hay muchos outliers por lo que se remueven para realizar el clustering.

mag5_7d_sinout <- mag5_7d[-which(mag5_7d$Aumento %in% outl_mag5), ]
mag5_14d_sinout <- mag5_14d[-which(mag5_14d$Aumento %in% outl_mag5), ]

mag6_7d_sinout <- mag6_7d[-which(mag6_7d$Aumento %in% outl_mag6), ]
mag6_14d_sinout <- mag6_14d[-which(mag6_14d$Aumento %in% outl_mag6), ]

mag7_7d_sinout <- mag7_7d[-which(mag7_7d$Aumento %in% outl_mag7), ]
mag7_14d_sinout <- mag7_14d[-which(mag7_14d$Aumento %in% outl_mag7), ]

par(mfrow=c(1, 3))

boxplot(mag5_7d_sinout$Aumento, mag5_14d_sinout$Aumento, main="Aumento para magnitud >= 5\ny 7 y 14 días sin outliers")
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

boxplot(mag6_7d_sinout$Aumento, mag6_14d_sinout$Aumento, main="Aumento para magnitud >= 6\ny 7 y 14 días sin outliers")
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

boxplot(mag7_7d_sinout$Aumento, mag7_14d_sinout$Aumento, main="Aumento para magnitud >= 7\ny 7 y 14 días sin outliers")
lines(x=c(0, 5), y=rep(0, 2), col="red", lwd=2)

Aumento ETA con respecto a Sismos sin outliers

Generación de los clusters con K-Means

Como se vio en la sección anterior, se aprecia que los datos están cargados hacia el positivo, es decir, hacia el aumento de los reportes de ETAs.

Dados los hallazgos encontrados en el preprocesamiento y la limpieza de datos, se procedió a realizar un clustering mediante el método de K-Means sobre los datos. El objetivo es averiguar si existen grupos de sismos que permitan predecir aumentos o disminuciones de reportes.

Se procede a trabajar para 7 y 14 días puesto que para estos valores se encontraron mejores resultados en la iteración anterior, y porque para los datos de 30 días la dispersión encontrada es muy amplia.

Generación de los datos:

generar_datos_cluster <- function(datos) {
  x <- data.frame(
    Magnitud = datos$Magnitud,
    Profundidad = datos$Profundidad,
    Region = datos$Region,
    Aumento = datos$Aumento
  )
  (scale(x))
}

x1 <- generar_datos_cluster(mag5_7d_sinout)
x2 <- generar_datos_cluster(mag5_14d_sinout)
x3 <- generar_datos_cluster(mag6_7d_sinout)
x4 <- generar_datos_cluster(mag6_14d_sinout)
x5 <- generar_datos_cluster(mag7_7d)
x6 <- generar_datos_cluster(mag7_14d)

Vemos la cantidad de clusters según el método del codo. Puesto que los datos de magnitud mayor o igual a 6 incluyen a los de 5, se procede con 5, a 7 días:

elbow_fun <- function (data) {
  wss <- 0
  for (i in 1:clust){
    wss[i] <- sum(kmeans(data, centers=i)$withinss)
  }
  wss
}

# Método del codo para 15 clusters, sobre el dataset para magnitud 6 y 7 días de intervalo
clust <- 15
wss1 <- elbow_fun(x1)
plot(1:clust, wss1, type="b", main="Codo para mag. 5, 7d.", xlab="Nº de clusters", ylab="wss")

Estrategia del codo para sismos y ETA

Se aprecia que el codo está en torno a los 5 clusters.

Se procede a realizar los clusters con K-Means. Para ello se ejecuta el siguiente script:

set.seed(3)
cant_clusters <- 5
ns <- 20

km.out1 <- kmeans(x1, cant_clusters, nstart=ns)
km.out2 <- kmeans(x2, cant_clusters, nstart=ns)
km.out3 <- kmeans(x3, cant_clusters, nstart=ns)
km.out4 <- kmeans(x4, cant_clusters, nstart=ns)
km.out5 <- kmeans(x5, cant_clusters, nstart=ns)
km.out6 <- kmeans(x6, cant_clusters, nstart=ns)

pairs(x1, col=km.out1$cluster, pch=8, main="Clustering para magnitud 5 y 7 días")
pairs(x2, col=km.out2$cluster, pch=8, main="Clustering para magnitud 5 y 14 días")
pairs(x3, col=km.out3$cluster, pch=8, main="Clustering para magnitud 6 y 7 días")
pairs(x4, col=km.out4$cluster, pch=8, main="Clustering para magnitud 6 y 14 días")
pairs(x5, col=km.out5$cluster, pch=8, main="Clustering para magnitud 7 y 7 días")
pairs(x6, col=km.out6$cluster, pch=8, main="Clustering para magnitud 7 y 14 días")

Se observa que para la región, no se encuentran clusters claros. La suposición hecha en el hito anterior de no utilizar la región para clasificar parece haber sido acertada. Se observa que para magnitud y profundidad, para las magnitudes 5 y 6 e intervalo de 7 y 14 días, se forman algunos grupos en que los ETAs aumentan y otros en que disminuyen.

Se podría decir que estos grupos sufren una especie de contrarrestación mutua, siendo el cluster con aumentos más claros el de color rojo generado para la magnitud 6, a 7 días de intervalo en la columna de Profundidad, y el color negro en la columna de Magnitud.

Clustering Ingresos Hospitalarions (IH)

Datos Sismológicos e Ingresos de pacientes a Hospitales

Los datos de IH son mucho mas numerosos que los de ETA por tanto veremos la posibilidad de analyzar solo los casos mas relevantes por Region. Bajo la asumpción de que los datos de ingresos de pacientes a centros hospitalarios incrementa cuando ocurren un evento sísmico, analizaremos las características principales de sismos con el numero de Ingresos a hospitales. El siguiente gráfico muestra la frecuencia de sismos en las regiones de Chile desde 2008 al 2018.

Note que solo consideramos el numero de IHs para definir el impacto de un sismo. Se omiten los demas datos que para este hito son considerados como relevantes. Corresponde a un trabajo futuro un análisis mas exhaustivo de estos datos y como afectan a la predicción.

library(ggplot2)

## load dataset of sismos
data <- read.csv("/Users/jhonc/Workspace/data-mining/back_r_sismos_chile/Datasets sismológicos/sismos_cleaned2.csv")
## query important sisms
sis_gt4 <- data[data$Magnitude >= 5 & data$year >= 2008 & data$year <= 2018, ]

ggplot(sis_gt4, aes(Magnitude)) + 
  # geom_col(position="dodge", stat="identity") 
  geom_bar() +
  facet_grid(year ~ Region_Number)

Matriz de histograma, Magnitud y Region

Vemos que los casos más notables de entre sismos y regiones se localizan entre las regiones 0 a 4, manteniéndose constante durante los años considerados. También podemos notar el caso especial de la región 13 el año 2010 y las regiones de 6 a la 9 también en el mismo año.

Consideremos para este análisis 3 casos en regiones 7ma en 2010, 1ra en 2014 y 4ta en 2015.

El siguiente codigo realiza un join de los datos, agregando la información de frequencia de IH en cada sismo. La agregación se realiza aplicando la suma.

library(dplyr)

egr_2010 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets egresoHospitales/cleaned2/egr_2010.csv", sep=",", header=TRUE)
egr_2014 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets egresoHospitales/cleaned2/egr_2014.csv", sep=",", header=TRUE)
egr_2015 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets egresoHospitales/cleaned2/egr_2015.csv", sep=",", header=TRUE)

ing_reg_10 <- egr_2010[egr_2010$REGION == 7,]
ing_reg_14 <- egr_2014[egr_2014$REGION == 1,]
ing_reg_15 <- egr_2015[egr_2015$REGION == 4,]

sis_10 <- with(sis_gt4, sis_gt4[Region_Number == 7 & year == 2010,])
sis_14 <- with(sis_gt4, sis_gt4[Region_Number == 1 & year == 2014,])
sis_15 <- with(sis_gt4, sis_gt4[Region_Number == 4 & year == 2015,])

## joined table
j_10 <- left_join(sis_10, ing_reg_10, by = c("number_days"))
j_14 <- left_join(sis_14, ing_reg_14, by = c("number_days"))
j_15 <- left_join(sis_15, ing_reg_15, by = c("number_days"))

col_list <- c('Depth', 'Magnitude', 'Latitude', 'Longitude', 'COMUNA','n')

j_10 <- j_10 %>%
  group_by(.dots=c("Latitude","Magnitude","COMUNA")) %>%
  add_tally()

j_14 <- j_14 %>%
  group_by(.dots=c("Latitude","Magnitude","COMUNA")) %>%
  add_tally()

j_15 <- j_15 %>%
  group_by(.dots=c("Latitude","Magnitude","COMUNA")) %>%
  add_tally()
# 

Numero de cluster optimos

En el siguiente gráfico mostramos los distintos valores del hiper-parámetro k en clustering para el grupo de datos seleccionados previamente. Note que para las regiones 7ma y 4ta estas tienen un descenso continuo, mientra que para 1ra region existe un salto en k=5. Sin embargo todos presentan mejora con el valor k=4, entonces definimos que k=4 es el número óptimo de clusters para todos los casos.

library(ggplot2, gridExtra)
library(cluster)
require(gridExtra)
library(cowplot)
################ kmeans 
set.seed(2)
clust = 15 # graficaremos hasta 15 clusters
elbow_fun <- function (data) {
  wss <- 0
  for (i in 1:clust){
    wss[i] <-
      sum(kmeans(data, centers=i)$withinss)
  }
  wss
}

col_list <- c('Latitude', 'Longitude', 'Depth', 'COMUNA', 'n')
# escalamos los datos, evitamos problemas a la hora de calcular el cluster
red_10 <- scale(j_10[col_list])
red_14 <- scale(j_14[!is.na(j_14$COMUNA),col_list])
red_15 <- scale(j_15[col_list])

wss1 <- elbow_fun(red_10)
wss2 <- elbow_fun(red_14)
wss3 <- elbow_fun(red_15)

par(mfrow=c(2,2))
p1 <- plot(1:clust, wss1, type="b", xlab="Numero de clusters (7ma Region) ", ylab="wss")
p2 <- plot(1:clust, wss2, type="b", xlab="Numero de clusters (1ra Region)", ylab="wss")
p3 <- plot(1:clust, wss3, type="b", xlab="Numero de clusters (4ta Region)", ylab="wss")

Tecnica del codo para datos sismos e IH

library(cluster)
# col_list <- c('Depth', 'Magnitude', 'Latitude', 'Longitude', 'COMUNA','n')
col_list <- c('Latitude', 'Longitude', 'Depth', 'Magnitude', 'n')

red_10 <- scale(j_10[col_list])
red_14 <- scale(j_14[!is.na(j_14$COMUNA),col_list])
red_15 <- scale(j_15[col_list])

cluster_10 <- kmeans(red_10, 4)
cluster_14 <- kmeans(red_14, 4)
cluster_15 <- kmeans(red_15, 4)

par(mfrow=c(2,2))
clusplot(red_10, cluster_10$cluster, color=TRUE, shade=FALSE, lines=1)
clusplot(red_14, cluster_14$cluster, color=TRUE, shade=FALSE, lines=1)
clusplot(red_15, cluster_15$cluster, color=TRUE, shade=FALSE, lines=1)
# plot_grid(p1, p2, p3)

Cluster de sismos e IH

Graficos de Correlación y clustering

El siguiente grafico muestra la matriz de cluster y correlación de la 7ma region en el año 2010, con k = 4. Note que para todos los casos, ‘n’ representa el numero de IH de un sismo en una fecha determinada.

Graficos de Correlación de atributos en clustering

El siguiente grafico muestra la matriz de cluster y correlación de la 7ma region en el año 2010, con k = 4. Note que para todos los casos, ‘n’ representa el numero de IH de un sismo en una fecha determinada.

d_10 <- j_10[col_list]
d_14 <- j_14[!is.na(j_14$COMUNA), col_list]
d_15 <- j_15[col_list]

fit <- kmeans(d_10, 8)
fit <- kmeans(d_14, 8)
fit <- kmeans(d_15, 8)

par(mfrow=c(3, 1))
pairs(d_10, col=fit$cluster)
pairs(d_14, col=fit$cluster)
pairs(d_15, col=fit$cluster)

Note que las correlaciones no son claras, la forma de los datos es densa. Sin embargo se pueden encontrar ciertas relaciones en la profundidad (Depth y la Longitud), un breve análisis muestra que esto es debido a la cercania de los sismos a las costas o dicho de una mejor manera al cinturon de fuego. El cinturon de fuego lleva como una capa superior a los continenetes que lo componen, entonces es mejor decir que a menor longitud mas profundiad por la forma de las capas de la tierra en esta zona del planeta. Esta relaciones es mejor representada en la region primera (Tarapaca) de Chile. No es el proposito de este estudio realizar mas conclusiones a estas relaciones de profundidad de un sismo.

Por otra parte, los grupos encontrados con respecto a la cantidad de IH forman grupos muy interesantes en la region cuarta. Los grupos se encuentran bien definidos en relación a la profundidad. Existe una gran cantidad ingresos hospitalarios a baja profundidad. Se puede concluir entonces que la predicción de incremento de incremento de IH con respecto a eventos sismicos no es trivial.

Clasificacion

Clasificadores de enfermedades de transmisión alimentaria

De acuerdo a lo observado en la sección anterior, se justifica realizar un clasificador, debido a que se logran apreciar algunos clusters relacionados con la profundidad y la magnitud del sismo y el aumento de los ingresos por enfermedades de transmisión alimentaria. Además, se observó que con sismos mayores o iguales a magnitud 6, los ingresos aumentaban y se mantenían en un valor positivo de diferencia respecto a la cantidad posterior al sismo y anterior a él. Con los magnitud 7 se observó un aumento de los ingresos post sismo siempre.

Para realizar una correcta clasificación, se sugirieron dos conjuntos de clases en base a la cantidad de «Aumento» planteado en las secciones anteriores.

2 clases: aumento o disminución

Clasificador binario. Debe predecir si un sismo produciría aumento en los ingresos o disminución, considerando el 0 como disminución.

Se clasifica de tal forma que si la columna «Aumento» es menor a 0, significa que disminuyó y, si no, aumentó.

4 clases: aumento/disminución alto/bajo.

Se establece una graduación a la cantidad de aumento, de tal forma que se distinga entre un aumento (o disminución) grande o bajo. Se justifica por los resultados del cluster.

Se considera bajo si el aumento o la disminución es menor a 12 (mayor a -12 en caso de ser disminución), mientras que alto en caso contrario.

En resumen, la clasificación queda como:

• Aumento < -12: disminución alta
• Aumento <= 0: disminución baja
• Aumento <= 12: aumento alto
• Caso contrario: aumento alto

La razón de clasificar con mayor cantidad de clases es despegarse de posibles variaciones no significativas, en este caso consideradas como variaciones menores a 12 ingresos pre y post sismo. De esta manera, se podría decidir de mejor forma si un sismo incide o no en la cantidad de ingresos post sismo. Además, el cluster arrojó al menos 4 tipos de sismos distintos, dos distinguibles entre aumento y disminución y dos que no tendrían relación con el aumento.

En la siguiente sección se explica por qué se escogió 12 como umbral.

Determinación del umbral para clasificador de ETAs

Para el nuevo clasificador, fue necesario determinar un umbral con tal de definir la disminución marcada, la disminución tenue, el aumento tenue y el aumento marcado.

Se cargan los datos con la función generar_datos, de la forma

mag4_7d <- generar_datos(sismos, etas, 4, 7, escribir=FALSE)

Luego se determinó el número para el umbral observando el comportamiento parecido a una distribución normal observado para los datos de magnitud 5

quantile(mag4_7d$Aumento, prob=seq(0, 1, length = 11))

Evaluación de la influencia de la región o no

Como hipótesis, se plantea que la región no debería influir en los resultados de la clasificación. Sin embargo, dado que en el norte se tienen más datos que en el sur de Chile (han habido más sismos y, por consiguiente, más sismos de magnitud mayor a 5), podría observarse una diferencia en las observaciones. Por esto, se realizaron experimentos de clasificación considerando las regiones y sin considerarlas.

Clasificador de árbol de decisión

Se decidió evaluar el desempeño de un árbol de decisión, con el que se intentaría predecir el impacto de los terremotos. Para poder observar las decisiones que va considerando el árbol de decisión, se buscó una forma de imprimir gráficamente el árbol, mostrando la iportancia de cada atributo.

Se eligió un árbol de decisión, debido a que las hojas pueden estar en diferentes grupos, sin importar la cantidad que el grupo tenga. Así, podría adaptarse a los distintos sismos.

Competencia de árboles de decisión

Cómo los resultados variaban dependiendo la profundidad con la que se formaban los árboles, se decidió hacerlos competir variando los parámetros profundidad máxima y criterion.

Competencia de clasificadores

Como el árbol no arrojó buenos resultados, como se muestra en la sección de resultados, se decidió probar otro tipo de clasificadores. Algunos, podrían ser más parecidos a los que se apreció en los clusters, como el KNN. Sin embargo, se decidió que no se limitaría a investigar el comportamiento de ese solamente, por lo que también se optó por usar un clasificador bayesiano y uno «dummy» que sería bueno para comparar el correcto entrenamiento de los clasificadores.

Resultados Clasificador para enfermedades de transmisión alimentaria

Resultados clasificador de árbol de decisión

Las siguientes librerías son necesarias para ejecutar los scripts:

import numpy as np
import pandas as pd

from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV

# Imports para el exportador de árbol
from sklearn.externals.six import StringIO
from IPython.display import Image
from sklearn.tree import export_graphviz
from sklearn import tree
import pydotplus

# Directorio donde están los datos para el clasificador
datos_dir= "Diferencia ETAs por magnitud/"

Clasificador binario

Este es el clasificador que se usó para las clases de «Aumento» o «Disminución».

def exportTreeImage(dtree, X, y, path):
    # Create dot data
    dot_data = tree.export_graphviz(dtree, out_file=None,
            filled=True,
            rounded=True, special_characters=True,
            feature_names=X.columns.tolist(),
            class_names=['Disminución', 'Aumento'])

    # Draw graph
    graph = pydotplus.graph_from_dot_data(dot_data)

    # Show graph
    Image(graph.create_png())

    # Create PNG
    graph.write_png(path)
    

def entrenarClasificador(prueba, path_X, path_y):
    print("########## Inicio de entrenamiento ###########")
    print(type(prueba))
    print("Leyendo archivo " + path_X)
    X= pd.read_csv(path_X);
    print("Success")
    print("Leyendo archivo " + path_y)
    y= pd.read_csv(path_y);
    print("Success")
    
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.33, random_state=9, stratify=y)

    clf= DecisionTreeClassifier(max_depth=5)
    clf.fit(X_train, y_train)

    y_pred= clf.predict(X_test)

    print("Accuracy en test set:", accuracy_score(y_test, y_pred))
    print(classification_report(y_test, y_pred))

    print("########### Fin de entrenamiento #############") 
    print("##############################################")

    return clf, X, y, X_train, y_train, X_test, y_test


def main():

    pruebas= ["5_7d", "6_7d", "7_7d", "5_14d", "5_7d_desh"]
    
    for p in pruebas:
    
        clf, X, y, X_train, y_train, X_test, y_test = entrenarClasificador(p,
                                                datos_dir + "X" + p + ".csv",
                                                datos_dir + "y" + p + ".csv")

        # Exportar árbol
        #exportTreeImage(clf, X, y, p + "binario.png")
main()

Clasificador de 4 clases

Este es el clasificador utilizado para las clases de Aumento/Disminución Alta/Baja.

def exportTreeImage(dtree, X, y, path):
    # Create dot data
    dot_data = tree.export_graphviz(dtree, out_file=None,
            filled=True,
            rounded=True, special_characters=True,
            feature_names=X.columns.tolist(),
            class_names=['Alta_Dism.', 'Baja_Dism.','Bajo_Aum.', 'Alto_Aum.'])

    # Draw graph
    graph = pydotplus.graph_from_dot_data(dot_data)

    # Show graph
    Image(graph.create_png())

    # Create PNG
    graph.write_png(path)
    

def entrenarClasificador(prueba, path_X, path_y):
    print("##############################################")
    print("########## Inicio de entrenamiento ###########")
    print(type(prueba))
    print("Leyendo archivo " + path_X)
    X= pd.read_csv(path_X);
    print("Success")
    print("Leyendo archivo " + path_y)
    y= pd.read_csv(path_y);
    print("Success")
    
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.33, random_state=9, stratify=y)

    clf= DecisionTreeClassifier(max_depth=5)
    clf.fit(X_train, y_train)

    y_pred= clf.predict(X_test)

    print("Accuracy en test set:", accuracy_score(y_test, y_pred))
    print(classification_report(y_test, y_pred))

    print("########### Fin de entrenamiento #############") 
    print("##############################################")

    return clf, X, y, X_train, y_train, X_test, y_test


def main():

    pruebas= ["5_7d", "5_14d", "6_7d", "6_14d", "7_7d", "7_14d"]

    for p in pruebas:
    
        print("##########################################")
        print("Clasificador con región")
        clf, X, y, X_train, y_train, X_test, y_test = entrenarClasificador(p,
                                                datos_dir + "X" + p + ".csv",
                                                datos_dir + "y" + p + "_cat.csv")

        # Exportar árbol
        #exportTreeImage(clf, X, y, p + ".png")

        print("##########################################")
        print("Clasificador sin región")
        clf, X, y, X_train, y_train, X_test, y_test = entrenarClasificador(p,
                                                datos_dir + "X" + p + "_sin_reg.csv",
                                                datos_dir + "y" + p + "_cat.csv")
        #exportTreeImage(clf, X, y, p + "sin_reg.png")
        
main()

Los resultados gráficos de los árboles presentados son los que se presentan a continuación. Se mostrarán solo algunos, pero los demás pueden verse en la carpeta «Diferencia ETAs por magnitud».

Aquí se pueden observar los resultdos distintos de árboles que se obtuvieron para los 14 días, con magnitudes de sismo mayores o iguales a 5. Se tienen muestras con los datos sin considerar las regiones, considerándolas y con 4 y 2 clases.

Árbol de decisión generado para sismos de magnitud mayor o igual a 5 con dos clases

Árbol de decisión generado para sismos de magnitud mayor o igual a 5 con 4 clases

Árbol de decisión generado para sismos de magnitud mayor o igual a 5 con 4 clases y sin considerar las regiones

Árbol de decisión generado para sismos de magnitud mayor o igual a 6 con 4 clases

Árbol de decisión generado para sismos de magnitud mayor o igual a 6 con 4 clases, sin considerar las regiones

Árbol de decisión generado para sismos de magnitud mayor o igual a 7 con 4 clases. Nota: en este caso, se imprimen mal los nombres de las clases, pues no se presentan todas. Alta disminución corresponde a Bajo aumento. Baja disminución corresponde a Alto aumento.

Árbol de decisión generado para sismos de magnitud mayor o igual a 7 con 4 clases. Nota: en este caso, se imprimen mal los nombres de las clases, pues no se presentan todas. Alta disminución corresponde a Bajo aumento. Baja disminución corresponde a Alto aumento.

Como norma general, se puede apreciar que las regiones se presentan muy arriba en el árbol. Es decir, se llevan la mayor cantidad del peso de las decisiones que este toma. Esto puede deberse a que la mayor cantidad de sismos de los que se tiene registro pueden alterar los resultados. El caso extremo fue analizar los sismo de magnitud mayor o igual a 7, en los cuales la región fue el único parámetro utilizado para clasificar los sismos. Por esta razón, se decidió corroborar la información con los mismos datos, pero sin considerar las regiones.

Sin las regiones, se obtuvieron árboles más pequeños en general, lo que habla también de una posible mejoría respecto a evitar sobreajustes.

A pesar de esto, como muestran los resultados de arriba, el árbol es bastante malo para encontrar las clases, en todos los casos. El accuracy y el recall son muy bajos para ser considerado un clasificador efectivo.

Se observa que la decisión de aumentar la cantidad de clases disminuyó el accuracy y el recall, por lo que no fue una medida efectiva.

El mejor de todos fue el clasificador que usó los datos de magnitud 6 para 7 días binario, como en el hito anterior.

Resultados de la competencia de clasificadores

Se tomaron magnitudes distintas para evaluar los clasificadores: 5, 6, y 7, de 14 días todos, que fueron los que tuvieron mejor capacidad predictiva según los resultados anteriores.

Además, se realizó competencia para aquellos binarios (con dos clases) y los con 4 clases. El código se encuentra en los anexos.

Binarios

Magnitud 5, 14 días

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.5057435931464518	0.5372547395781313	0.5778510470900172	0.5358593799236964
Recall promedio	0.5035205992509363	0.5334831460674156	0.5755430711610486	0.5362921348314607
F1-score promedio	0.5032128058844465	0.5334884367150398	0.565797022911816	0.5345539655623789

Con 4 clases

Magnitud 5, 14 días

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.2543251271511943	0.28198277519152576	0.28500778188199705	0.29402535668893387
Recall promedio	0.25168539325842704	0.2779026217228464	0.3156554307116105	0.2847191011235955
F1-score promedio	0.2513672342633679	0.277833332404425	0.2559771146372389	0.28464672910000527

Magnitud 6, 14 días

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.31734917970917964	0.29914288933288935	0.2681247261059187	0.220177534337008
Recall promedio	0.2912	0.2716	0.28519999999999995	0.18959999999999996
F1-score promedio	0.28493039675184	0.26893306454700533	0.23373182632714543	0.1826686267078818

Magnitud 7, 14 días

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.5138888888888888	0.8622222222222221	0.6283333333333334	0.4461111111111111
Recall promedio	0.43000000000000005	0.8433333333333333	0.6366666666666667	0.5033333333333333
F1-score promedio	0.41966666666666663	0.8326666666666666	0.6016666666666667	0.43833333333333324

Se observa que, en general, todos los clasificadores se asemejan al clasificador Dummy en los resultados. Se concluye, entonces que los clasificadores no fueron una buena herramienta, probablemente debido a la falta de datos (en general fueron pocos, por la falta de sismos fuertes) o a distintos comportamientos entre regiones. En este sentido, para el futuro podrían aislarse las regiones y probar por cada una de ellas, para que el clasificador se «adapte» a su región en particular.

Clasificador para Ingresos Hospitalarios

El desarrollo del experimento de clasificación en IH se realizo de manera analoga a la de ETA. A continuación mostraremos el desarrollo de este expermiento.

Distribución de los incrementos/decrementos de Ingresos Hospitalarios

La distribusion de diferencia entre incremento y decremento muestra una distribusion quasi-normal. Donde las colas se encuentran estiradas e inclinadas a los decrementos. Esta diferencia no es mas que una simple suma entre cambios del dia en el que se produjo un sismo \(d - delta\) + \(d + delta\). Recordando siempre que se toman los sismos >= 5 y tomando como delta = {3, 7, 14}. Recuerde que cuando analisamos los sismos de >= 6 o >=7 son solo un sub-conjunto de este.

d1 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets sismológicos/ing_data/ing_sis5_d3.csv")
d2 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets sismológicos/ing_data/ing_sis5_d7.csv")
d3 <- read.csv("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets sismológicos/ing_data/ing_sis5_d14.csv")

par(mfrow=c(2,2))
plot(density(d1$ING_DIFF),  main="Delta=3")
rug(jitter(d1$ING_DIFF))

plot(density(d2$ING_DIFF),  main="Delta=7")
rug(jitter(d2$ING_DIFF))

plot(density(d3$ING_DIFF),  main="Delta=14")
rug(jitter(d3$ING_DIFF))

boxplot(d1$ING_DIFF,d2$ING_DIFF,d3$ING_DIFF, names=c("3","7","14"))

Densidad para sismos >= 5 y diferencia entre incremento y decremento con delta = {3,7,14}

Clasificador binario

Clasificamos los datos binarios de los sismos, como incremento o decremento. Utilizando el mismo análisis que en ETA. Obtenemos los datos utilizando un simple algoritmo utilizando la diferencia ya mencionada previamente.

La distribusion mostrada anteriormente nos dice que los datos no se encuentran del todo balanceados. Por tanto haremos un balance de datos utilizando la estrategia oversampling o upsampling, donde incrementamos el numero de la clase que se encuentran con menor tamaño.

Resultados de clasificadores con magnitud >= 5 y delta = 3

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.4045992532391457	0.43423935482523673	0.4312893109382294	0.43763807844250396
Recall promedio	0.48902361685024626	0.38223832381772865	0.8143813666109686	0.5045940698252863
F1-score promedio	0.4421608273227311	0.4059665700288416	0.5600709255119732	0.4681348468758294

Resultados de clasificadores con magnitud >= 5 y delta = 7

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.3872961393667831	0.4193564490530995	0.4110486648046686	0.41380027174033207
Recall promedio	0.4954397469934727	0.3648461502302095	0.9629178373452127	0.48359961538226237
F1-score promedio	0.43426467191138485	0.3894902181009887	0.5757260418739462	0.44538193028386097

Resultados de clasificadores con magnitud >= 5 y delta = 14

Clasificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.4015495834059787	0.42620754088348145	0.41532882054787024	0.43428140414898986
Recall promedio	0.5050258066525183	0.3782505652987394	0.949464422546	0.48749837079762254
F1-score promedio	0.44706266563950176	0.3998933267205183	0.5773383527694784	0.45843017970970346

Se puede apreciar que el los clasificadores no son tienen un buen resultado en precision superando por poco a DummyClassifier. Dentro de los clasificadores tambien se puede apreciar que el clasificador NaiveBayes es el que tiene mejores resultados en recall en promedio. Este resultado es interesante puesto que este problema predice información sensible. Es decir que nos favorece tener predicción sobre falsos positivos y no en falsos negativos. Esta sentencia puede ser facilmente respaldada proponiendo un ejemplo hipotetico, si un evento sismico sucede en un 18 de Septiembre en la noche de la mayoria de la pobracio chilena se encuentra, ya en un evento o en sus casas, es mejor que el predictor haga mension que el incremento de ingresos hospitalarios, asi poder estar listos ante cualquier contingencia.

Classificación Multi-clase

El clasificador multi-clase fue configurado analogamente a el análisis de ETA. Primero calculamos los cuantiles de los datos que tenemos.

quantile(d1$ING_DIFF, prob=seq(0, 1, length = 11))
summary(d1$ING_DIFF)

La diferencia en ingresos se encuentra inclinada negativamente. En este caso la información mas sensible se encuentra en los incrementos, que es lo que desamos clasificar y predecir en un futuro. Entonces tomamos los el cuantil positivo (3rd cuantil) y utilizamos un umbral de valor = 75.

En resumen, la clasificación queda como:

• Aumento < -75: disminución alta
• Aumento <= 0: disminución baja
• Aumento <= 75: aumento bajo
• Caso contrario: aumento alto

Resultados de clasificadores con magnitud >= 5 y delta = 3

Classificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.2507327576140339	0.3876901299003524	0.5058967666134677	0.3543702834338281
Recall promedio	0.25038281649044	0.39663779445334424	0.44343383481408	0.36345213837534623
F1-score promedio	0.2439288784869155	0.38944852613051834	0.3789695690938443	0.3546972245298021

Resultados de clasificadores con magnitud >= 5 y delta = 7

Classificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.24512935121337415	0.3583267783121016	0.4804950794067688	0.33050904858478936
Recall promedio	0.24634647647219773	0.35713542142458576	0.35688271114936576	0.32835525998190407
F1-score promedio	0.24230364232042562	0.3554661615286645	0.2992047260360324	0.3274396517625613

Resultados de clasificadores con magnitud >= 5 y delta = 14

Classificador	Base Dummy	Decision Tree	Gaussian Naive Bayes	KNN
Precision promedio	0.2563906674778534	0.3624337639804715	0.39474347014438027	0.3198595637440131
Recall promedio	0.25317694415974595	0.3686265697204939	0.39462598654478326	0.3293566817345862
F1-score promedio	0.24580116039825536	0.3629313448757567	0.3230796805553361	0.3167516076350386

Aqui los resultados resultan satisfactorios. Si bien anteriormente apreciamos que NaiveBayes tiene buenos resultados en recall en promedio, en este caso no podemos decir lo mismo.

Analysis de Clasificación utilizando Decision-Tree

Los arboles de decision nos pueden revelar información relacionada al creterio de clasificación cuando las relaciones no son muy claras. En el siguiente análisis vemos que los criterios son los esperados a que un sismo tenga un gran impacto o no. Es decir la magnitud y profundidad.

La imagen de abajo representa a un árbol de profundidad 1 utilizando el algoritmo DecisionTree. Cuando el árbol tiene una profundidad 1 el criterio de clasificación es la profundidad, lo cual representa la predominancia de este atributo para la predicción, igual como se mostró en el análisis de los datos.

Resultado del algoritmo DesicionTree con profundidad 1

Resultado del algoritmo DesicionTree con profundidad 3

Conclusiones

Estudiar la relación de sismos con otros eventos ayuda a la toma de decisiones que mitigan los efectos negativos. Gracias a la gran cantidad de información recolectada de sismos en los ultimos años es posible su mejor estudio. Sin embargo es complicado analyzar su relación con otros eventos por falta de información satelite publica. Dentro de este análisis estudio las posibles correlaciones entre sismos y eventos humanos, como Enfermedades de Transmision alimentaria (ETA) e Ingresos Hospitalarios (IH). Este estudio busca explotar estas correlaciones de estos para poder crear un predictor ayude a detectar incremento en ETA e IH. Llegamos a ver una leve correlación entre estos datos, gastamos varios recursos de data-mining para explorar la data. Tuvimos mucho cuidado manipulando la data y no introducir un error en el análisis.

Creamos un algorimo para asignar regiones dentro de Chile a los eventos sismicos, el cual no se encuentra publico en ninguna parte. Aprovechamos la geografia chilena para generar delimitaciones que reflejen las delimitaciones geograficas de las regiones. Filtramos información que nos fue util para el análisis. Consecuentemente utilizamos el tiempo y region chilena para poder unir los datos satelite, ETA e IH. Donde utilizamos tecnicas de agregación para poder medir y unir la información satellite a los sismos. Este estudio se enfoco al volumen de estos eventos, tomadolos como echos atomicos. Es decir, no se considera otra información que venia junto a estos eventos humanos, ETA e IH. El corto tiempo de esta información no permitio ser mas estricto en este aspecto, es de un trabajo futuro poder ser mas explicitos en esta informacion.

Durante el análisis de correlación utilizando clustering y correlación lineal vimos cierta relación entre la profundidad, magnitud y el incremento de ETA e IH, como tambien la region. Esto esta directamente relacionado al impacto de un sismo que usualmente se mide con la escala de Mercalli. Data que mida los sismos con la escala de Mercalli no pudo ser encontrada por ningun lado. Recordar que esta solo es encontrada cuando existen eventos que generen impacto negativo economico o humano, ademas de ser una escala subjetiva. La relación de eventos humanos como ETA e IH son relacionados con el impacto de un sismo y este tambien depende a que tan bien la localidad del evento se encuentre preparado. Entonces la region donde ocurre un sismo tambien influye en los eventos humanos.

Por ultimo creamos dos tareas de clasificación , binaria y multi-clase. Los cuales buscan predecir incrementos de eventos humanos como ETA e IH. Los resultados en general no son buenos, comparados con un clasificador de juguete como DummyClassifier. Interesantemente encontramos una sorpresa en el clasificador binario, el clasificador Naive Bayes tuvo un buen score en recall en promedio, llegando a tener un recall de 0.96%. La medida de recall es utilizada cuando se una tarea que incluye predicción que es sensible. Es decir, que predice mas veces en predecir una falsa alarma. En este caso resulta ser considerado como positivo ya que ayuda a tener planes de contingencia ante un evento sismico.

Anexo

Un simple análisis de sismos y precipitaciones

Se busca encontrar una relación entre precipitaciones y sismos. En este análisis se toma encuenta los datos de temporadas de año, son temporadas percivibles en las cuales podemos facilmente decernir si estos mantienen una relacion. Ademas buscamos relación de sismos y los meses para tener un resultado mas cercano a fechas concretas.

Sismos por estaciones del año

La siguiente es una gráfica de sismos a través de las estaciones del año tomando en cuenta datos desde el año 1960 hasta el 2019. La grafica muestra que existe una aparente correlación entre las temporadas del año y la cantidad de terremotos.

Relación lluvia y sismos

Relación entre precipitaciones y sismos por mes.

Los datos fluviales fueron extraídos de la página web de CR2. Fueron tomados desde enero de 1900 hasta febrero de 2018, en 874 estaciones de todo Chile. Existen 2 variantes: la resolución temporal mensual o diaria. Las precipitaciones están medidas en milímetros. A continuación se muestra las caracteristicas del dataset, y el tratamiento y limpieza de los datos.

Las dimensiones del dataset corresponden a 874 filas y 1431 columnas, correspondientes a los datos obtenidos en todos los años desde 1900 a 2018, por cada estación. Además, presentan datos como la latitud y longitud de la estación, la altura a la que se encuentra sobre el nivel del mar y la cuenca a la que pertenecen.

Este dataset se dejó de usar para el Hito 3, porque en las iteraciones anteriores no se encontraron correlaciones entre los sismos y las lluvias, por lo que no resultó relevante utilizar técnicas de minería de datos para revisar estos dataset.

# Importar dataset de lluvias mensuales
prAmon <- t(read.csv("https://anakena.dcc.uchile.cl/~cllull/IntroMineriaDatos/DataSets/cr2_prAmon_2018/cr2_prAmon_2018.txt"))

colnames(prAmon) <- as.character(unlist(prAmon[1, ])) # Le pone nombre a las columnas
prAmon <- prAmon[-1, ] # Extrae los datos

prAmon_na <- prAmon # Copia de prAmon
prAmon_na[prAmon == -9999] <- NA # Todas las celdas con -9999 a NA

Se deben declarar como NA los datos con -9999 que representan datos faltantes, para que R los trabaje de buena manera.

Luego se aplica el promedio de todos los años desde 1900 a 2018, evitando los NA.

library(dplyr)
data <- prAmon_na[, 15:1431]  # Los datos de lluvia, sin descripción
data <- as.data.frame(data)  # De matrix a data.frame
data[,] <- apply(data[,], 2, function(x) as.numeric(x)) # Valores de character a numeric
prAmonMean <- prAmon_na[, 1:14]  # Tabla para poner el promedio
prAmonMean <- as.data.frame(prAmonMean)  # De matrix a data.frame
mean <- rowMeans(data, na.rm = TRUE)  # Calcula el promedio
prAmonMean$mean <- mean  # Pone el promedio en la tabla

Se separaron los datos por región, debido a que Chile tiene muchos climas diferentes.

prAmonMean$latitud <- as.numeric(as.character(prAmonMean$latitud)) # Se debe hacer numérico
prAmonMean_arica <- filter(prAmonMean, latitud <= -17.46 & latitud >= -19.07)
prAmonMean_magallanes <- filter(prAmonMean, latitud < -49.10 & latitud >= -56)

Lo más ilustrativo que se puede obtener es determinar las precipitaciones en cada mes del año. De la siguiente forma se extraen las precipitaciones mensuales para la región de Arica:

Precipitaciones en Arica y Parinacota

A continuación una grafico de los sismos por mes vemos que estos no tiene ninguna relación con el grafico de Promedio de precipitaciones por mes. Aunque el grafico es algo mas simple ese muestra no concordancia

# load libraries
library(tidyverse)
library(lubridate)

ggplot(season_data, mapping = aes(x=month)) +
  geom_bar() + 
  xlab("Mes") + ylab("Promedio precipitaciones") + coord_flip(expand = TRUE)

Precipitaciones en Arica y Parinacota

Codigo Fuente

# Agregarle las fechas a las ETAs
agregarFechasEtas <- function(dataframe) {
  fechas <- as.Date(dataframe$Fecha.de.Ingestión, "%d-%m-%Y")
  fechas <- data.frame(Year=as.numeric(format(fechas, format="%Y")),
                       Month=as.numeric(format(fechas, format="%m")),
                       Day=as.numeric(format(fechas, format="%d")))
  fecha_num <- data.frame(num_date=fechas$Year +
                            (fechas$Month-1)/12 +
                            (fechas$Day-1)/24/12)
  fechas <- cbind(fechas, fecha_num)
  cbind(dataframe, fechas)
}

#################################
# Competencia de clasificadores #
#################################
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import f1_score, recall_score, precision_score


def run_classifier(clf, X, y, num_tests=100):
    metrics = {'f1-score': [], 'precision': [], 'recall': []}
    
    for _ in range(num_tests):
        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.30)
        
        clf.fit(X_train, y_train)
        
        predictions = clf.predict(X_test)
        
        metrics['y_pred'] = predictions
        metrics['y_prob'] = clf.predict_proba(X_test)[:,1]
        metrics['f1-score'].append(f1_score(y_test, predictions, average='weighted')) 
        metrics['recall'].append(recall_score(y_test, predictions, average='weighted'))
        metrics['precision'].append(precision_score(y_test, predictions, average='weighted'))
    
    return metrics

import pandas as pd

from sklearn.dummy import DummyClassifier
from sklearn.svm import SVC  # support vector machine classifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.naive_bayes import GaussianNB  # naive bayes
from sklearn.neighbors import KNeighborsClassifier


path_X = datos_dir + "X7_14d.csv"
path_y = datos_dir + "y7_14d_cat.csv"

print("Leyendo archivo " + path_X)
X= pd.read_csv(path_X);
print("Success")
print("Leyendo archivo " + path_y)
y= pd.read_csv(path_y);
print("Success")

c0 = ("Base Dummy", DummyClassifier(strategy='stratified'))
c1 = ("Decision Tree", DecisionTreeClassifier())
c2 = ("Gaussian Naive Bayes", GaussianNB())
c3 = ("KNN", KNeighborsClassifier(n_neighbors=5))

classifiers = [c0,c1, c2, c3]

results = {}
for name, clf in classifiers:
    metrics = run_classifier(clf, X, y)   # hay que implementarla en el bloque anterior.
    results[name] = metrics
    print("----------------")
    print("Resultados para clasificador: ",name) 
    print("Precision promedio:",np.array(metrics['precision']).mean())
    print("Recall promedio:",np.array(metrics['recall']).mean())
    print("F1-score promedio:",np.array(metrics['f1-score']).mean())
    print("----------------\n\n")

Codigo de asignación de Regiones para sismos

separateRegions <- function(DataFrame){
  DataFrame$Region_Arica_Parinacota <- 0
  DataFrame$Region_Tarapaca <- 0
  DataFrame$Region_Antofagasta <- 0
  DataFrame$Region_Atacama <- 0
  DataFrame$Region_Coquimbo <- 0
  DataFrame$Region_Valparaiso <- 0
  DataFrame$Region_Metropolitana_Santiago <- 0
  DataFrame$Region_Libertador_General_OHiggins <- 0
  DataFrame$Region_Maule <- 0
  DataFrame$Region_Nuble <- 0
  DataFrame$Region_Biobio <- 0
  DataFrame$Region_Araucania <- 0
  DataFrame$Region_Rios <- 0
  DataFrame$Region_Lagos <- 0
  DataFrame$Region_Aysen_General_Carlos_Ibanez <- 0
  DataFrame$Region_Magallanes_Antartica <- 0
  for (i in 1:nrow(DataFrame)) {
    if (DataFrame[i, "Region_Number"] == 1) {
      DataFrame[i, "Region_Tarapaca"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 2) {
      DataFrame[i, "Region_Antofagasta"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 3) {
      DataFrame[i, "Region_Atacama"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 4) {
      DataFrame[i, "Region_Coquimbo"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 5) {
      DataFrame[i, "Region_Valparaiso"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 6) {
      DataFrame[i, "Region_Libertador_General_OHiggins"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 7) {
      DataFrame[i, "Region_Maule"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 8) {
      DataFrame[i, "Region_Biobio"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] ==9) {
      DataFrame[i, "Region_Araucania"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 10) {
      DataFrame[i, "Region_Lagos"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 11) {
      DataFrame[i, "Region_Aysen_General_Carlos_Ibanez"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 12) {
      DataFrame[i, "Region_Magallanes_Antartica"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 13) {
      DataFrame[i, "Region_Metropolitana_Santiago"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 14) {
      DataFrame[i, "Region_Rios"] <- 1
    }
    else if (DataFrame[i, "Region_Number"] == 15) {
      DataFrame[i, "Region_Arica_Parinacota"] <- 1
    }
    else {
      DataFrame[i, "Region_Nuble"] <- 1
    }
  }
  
  return (data.frame(DataFrame))
}

Codigo de agregación de datos para Ingresos hospitalarios

calculate_prev_after <- function(reg, ing, k) {
  min <- reg[['number_days']] - k
  max <- reg[['number_days']] + k
  
  ing_data <-  with(ing, ing[REGION == reg[['Region_Number']], c('number_days', 'DIAG1', 'REGION','COMUNA')])
  ings_bef <-  with(ing_data, ing_data[min <= number_days & number_days < reg[['number_days']], ]) 
  ings_aft <-  with(ing_data, ing_data[reg[['number_days']] < number_days & number_days <= max, ]) 
  
  ings_bef <- nrow(ings_bef)
  ings_aft <- nrow(ings_aft)

  
  return(data.frame(ING_BEFORE = ings_bef, ING_AFTER = ings_aft, ING_DIFF= (ings_aft - ings_bef)))
}

add_columns <- function(dataset, window) {
  classes <- list()
  first <- 1
  for ( i in 1:nrow(dataset)) {
    # print (i)
    register <- dataset[i,]
    ds <- get_dataset(register[['year']])
    if (first == 1 ) {
      ing_data <- calculate_prev_after(register, ds, window)
      first <- 0
    } else {
      row <- calculate_prev_after(register, ds, window)
      ing_data <- rbind(ing_data, row)
    }
    
  }
  # classes
  ing_data
}

create_ing_data <- function(min_mag, window, data, file_path) {
  sis_gt4 <- data[data$Magnitude >= min_mag & data$year >= 2008 & data$year <= 2018, ]
  ing_data <- add_columns(sis_gt4, window)
  
  sis_gt4$ING_BEFORE <- ing_data$ING_BEFORE
  sis_gt4$ING_AFTER <- ing_data$ING_AFTER
  sis_gt4$ING_DIFF <- ing_data$ING_DIFF
  write.csv(sis_gt4, file = file_path)
  sis_gt4
}

codigo para generar data para clasificación binaria

for_binary_classification <- function (src_file, target_file) {
  d <- read.csv(src_file)
  d <- separateRegions(d)
  d <- transform(d, class=ifelse(ING_DIFF > 0, 1, 0))
  write.csv(d, file = target_file)
  d
}

codigo para generar data para clasificación con 4 clases

for_four_class_classification <- function (src_file, target_file) {
  d <- read.csv(src_file)
  d <- separateRegions(d)
  delta <- 75
  d <- transform(d, class=ifelse(ING_DIFF > 0, ifelse(ING_DIFF>delta, 3, 2), ifelse(ING_DIFF < -delta, 0, 1) ))
  write.csv(d, file = target_file)
  d
}

for_four_class_classification("/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets sismológicos/ing_data/ing_sis7_d3.csv",
      "/Users/jhonc/Workspace/data-mining/r_sismos_chile/Datasets sismológicos/four_class/test.csv")